module TriangleElementManipulation
    use MatrixOperation
    use GaussianIntegration
    use VectorOperations
    implicit none

    integer, parameter :: nNodes = 3
    integer, parameter :: probDim = 2 !< Problem Dimension | Plane Problem

    integer, parameter :: MatrixSize = nNodes * (probDim + 1) !< Heat Fluxes (each dimension) and Temperature problem
	integer, parameter :: XAxis = 1
	integer, parameter :: YAxis = 2

	integer, parameter :: HeatFluxX_Initial = 1 !< Keeping track of Degrees of Freedom index meaning
	integer, parameter :: HeatFluxX_Final   = HeatFluxX_Initial + nNodes - 1
	integer, parameter :: HeatFluxY_Initial = HeatFluxX_Final + 1
	integer, parameter :: HeatFluxY_Final = HeatFluxY_Initial + nNodes - 1
	integer, parameter :: Temperature_Initial = HeatFluxY_Final + 1
	integer, parameter :: Temperature_Final = Temperature_Initial + nNodes -1

	integer, parameter :: Dof_Values = 2
	integer, parameter :: Dof_Try = 1 !< Try values for the element DoF
	integer, parameter :: Dof_Commit = 2 !< Commit values for the element DoF

    !> \brief Triangular Element
    !! \warn The element is not aware of its Neighbours. Is up to the grid decide which is neighbour of which
    !! The Solution of the Discontinuous Finite Element formulation is by Sucessive Substitution, which is equivalent to
    !! iteratively solve a block-diagonal global matrix. This kind of solution is well-suited to parallel environments
    !! with fast interprocess communication systems. Due to this choice of solution procedures the local values of the
    !! Degrees of Freedom must be separated into two categories, the trying value and the committed value.
    type TriangleElement
!// \todo (ronaldo#1#): Modificar Matrix, Dof, OuterVector, NeighbourReNumbering para allocatable e ajustá-los em SetNodesTriangularElement
        real(kind=8), dimension(:,:), allocatable :: Matrix
        integer, dimension(nNodes)                :: Nodes
        real(kind=8), dimension(:,:)              :: DoF !< Degres-Of-Freedom are comprised for the Try and Commit Values
        real(kind=8), dimension(:,:)              :: OuterVector !< Outer vector oposed to the vertex
        integer, dimension(nNodes, nNodes)        :: NeighbourReNumbering
        real(kind=8) 							  :: S !< Auxiliar Space variable used to make calculations faster
    contains
        procedure, public :: SolutionErrorRange
        procedure, public :: SetNodesTriangularElement
        procedure, public :: SetMatrixTriangleElement
        procedure, public :: SetVectorTriangleElement
        procedure, public :: SetAllNeighboursNumbering

    end type

contains


    !> \brief Returns the Error range estimation for the current Element
    !!
    !! \param this - the element
    !! \return the maximum error over the components of the Degrees-Of-Freedom vector
    !!
	pure function SolutionErrorRange( this ) result (v)
		implicit none
		type(TriangleElement), intent(in) :: this
		real(kind=8) :: v
		v = norm_infty(this%DoF(:,Dof_Commit)-this%DoF(:,Dof_Try))
	end function

    !> \brief Set Nodes Indexes of a Triangular Element
    !! These are the nodes of the grid
    !! \param this - the element itself
    !! \param Nodes - the Indexes of the Nodes on the Grid
    !!
    subroutine SetNodesTriangularElement( this, Node1, Node2, Node3)
        implicit none
        type(TriangleElement), intent(out) :: this
        integer, intent(in) :: Node1
        integer, intent(in) :: Node2
        integer, intent(in) :: Node3


		integer :: statReturn

        this%Nodes(1) = Nodes1
        this%Nodes(2) = Nodes2
        this%Nodes(3) = Nodes3

		allocate( this%Matrix(MatrixSize, MatrixSize+1), stat=statReturn )
		if (statReturn /= 0) then
			print *, 'Unable to allocate memory for Matrix on Triangle Element.'
			return
		end if

        this%Matrix = 0.0d0

        allocate( this%DoF(MatrixSize, Dof_Values), stat=statReturn )
        if (statReturn /= 0) then
        	print *,  'Unable to allocate memory for Degrees of Freedom on Triangle Element.'
        	return
        end if

        this%DoF = 0.0d0

		allocate( this%NeighbourReNumbering(nNodes, nNodes), stat=statReturn )
		if (statReturn /= 0) then
			print *, 'Unable to allocate memory for Neighbour Renumbering Matrix on Triangle Element.'
			return
		end if

		this%NeighbourReNumbering = 0

        allocate( this%OuterVector(nNodes, probDim), stat=statReturn )
        if (statReturn /= 0) then
        	print *, 'Unable to allocate memory for Outer Vector on Triangle Element.'
        	return
        end if

    end subroutine





    !> \brief Groups the individual constructions of the several components of the Linear Systems Matrix
    !!
    !! \param this - The Triangular Element itself
    !! \param
	!!
	subroutine SetMatrixTriangleElement( this, X, ThermCond, Stabilization )
        implicit none
        type(TriangleElement), intent(inout) :: this
        real(kind=8), dimension(nNodes, probDim), intent(in) :: X
        interface
			pure function ThermCond(T) result (v) !< Thermal Conductivity Function
				implicit none
				real(kind=8), intent(in) :: T
				real(kind=8) :: v
			end function
		end interface
		real(kind=8), intent(in) :: Stabilization

		this % Matrix(1:MatrixSize,1:MatrixSize) = 0.0d0

        call SetMatricesE( this%Matrix(1:MatrixSize,1:MatrixSize), X )
        call SetMatricesH_l( this, this%Matrix(1:MatrixSize,1:MatrixSize), X )
        call SetMatricesJ_l( this, this%Matrix(1:MatrixSize,1:MatrixSize), X )

		call SetMatricesE_lr_i( this, this%Matrix(1:MatrixSize,1:MatrixSize), X, ThermCond, Stabilizaiton )
		call SetMatricesH_l_i( this, this%Matrix(1:MatrixSize,1:MatrixSize), X, ThermCond, Stabilizaiton )
		call SetMatricesJ_l_i( this, this%Matrix(1:MatrixSize,1:MatrixSize), X, Stabilizaiton )
		call SetMatrixG_T_i( this%Matrix(1:MatrixSize,1:MatrixSize), X, Stabilizaiton )
    end subroutine


	subroutine SetVectorTriangleElement( this, X, Neighbours, ThermCond, Stabilization )
        implicit none
        type(TriangleElement), intent(inout) :: this
        real(kind=8), dimension(nNodes, probDim), intent(in) :: X
        type(TriangleElement), dimension(nNodes), intent(in) :: Neighbours
        interface
			pure function ThermCond(T) result (v) !< Thermal Conductivity Function
				implicit none
				real(kind=8), intent(in) :: T
				real(kind=8) :: v
			end function
		end interface
		real(kind=8), intent(in) :: Stabilization

		call SetMatricesE_lr_B_i(this, this%Matrix(:,MatrixSize+1), Neighbours, ThermCond, Stabilization )
		call SetMatricesH_l_B_i(this, this%Matrix(:,MatrixSize+1), Neighbours, Stabilization, ThermCond )
		call SetMatricesJ_l_B_i(this, this%Matrix(:,MatrixSize+1), Neighbours, Stabilization )
	end subroutine

	subroutine SetAllNeighboursNumbering(this, Neighbours)
		implicit none
		type(TriangularElement), intent(inout) :: this
		type(TriangularElement), dimension(nNodes), intent(in) :: Neighbours

		integer, intent(in) :: i

		do i = 1, nNodes
			call SetNeighbourNumbering(this, i, Neighbours(i))
		end do
	end subroutine


	subroutine SetNeighbourNumbering(this, n, Neighbour)
		implicit none
		type(TriangularElement), intent(inout) :: this
		integer, intent(in) :: n
		type(TriangularElement), intent(in) :: Neighbour

        integer :: k, f
        do k = 1, nNodes
            f = findNode(k)
            if (f == 0) then
                this%NeighbourReNumbering(n,k) = k
            else
                this%NeighbourReNumbering(n,findNode(this, k, Neighbour)) = k
            end if
        end do
    end subroutine

    pure function findNode(this, n, Neighbour) result (NbIdx)
        implicit none
		type(TriangularElement), intent(inout) :: this
        integer, intent(in) :: n
		type(TriangularElement), intent(in) :: Neighbour

        integer :: NbIdx
        integer :: k
        do k = 1, nNodes
            if (Neighbour%nodes(k) == this%Nodes(n)) then
                NbIdx = k
                return
            end if
        end do
        NbIdx = 0
    end function


	pure function VectorRenumber( this, i ) result (Indexes)
		implicit none
		type(TriangularElement), intent(in) :: this
		integer, intent(in) :: i
		integer, dimension(MatrixSize) :: Indexes

		integer :: k
		do k = 1, nNodes
			Indexes(HeatFluxX_Initial+k-1) = HeatFluxX_Initial+this%NeighbourReNumbering(i, k)-1
			Indexes(HeatFluxY_Initial+k-1) = HeatFluxY_Initial+this%NeighbourReNumbering(i, k)-1
			Indexes(Temperature_Initial+k-1) = Temperature_Initial+this%NeighbourReNumbering(i, k)-1
		end do
	end function

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Calculations Zone Delimiter!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


	!> \brief Calculate the rotation of nodes over a Triangle
	!!
	!! \param n - the edge oposed to vertex n which is going to be calculated the normal vector
	!! \return v - 3-d integer vector with the local indexes of the Nodes
	!!
	pure function rotation(n) result (v)
		implicit none
		integer, intent(in) :: n
		integer, dimension(3) :: v
		select case(n)
			case(1)
				v(1) = 1
				v(2) = 2
				v(3) = 3
			case(2)
				v(1) = 2
				v(2) = 3
				v(3) = 1
			case(3)
				v(1) = 3
				v(2) = 1
				v(3) = 2
			case default
				v(1) = 1
				v(2) = 2
				v(3) = 3
		end select
	end function


    !> \brief Calculate the outer normal vector oposed to de vertex i
    !!	the calculations are preformed using a extension of the problem dimension to generate a orthogonal vector by cross vector product.
    !! \param this - Triangular Element
    !! \param Coordinates - The Coordinates Position of the Nodes
    !!
    subroutine CalculateOuterVector(this, Coordinates)
        implicit none
        type(TriangleElement), intent(out) :: this
        real(kind=8), dimension(nNodes,probDim) :: Coordinates

		real(kind=8), dimension(2) :: Xm
		real(kind=8), dimension(3) :: V1, V2

        integer :: n
		integer, dimension(3) :: p

        do n=1,3
			p = rotation(n)

			!< Xvirt = Xm(extended); Xm(3) = 1.0d0
			!< The projection of Xvirt on the plane is already over the the edge defined by Xp(2) and Xp(3)
			Xm = (Coordinates(p(2),:) + Coordinates(p(3),:))/2.0d0
			V1(1:2) = Xm-Coordinates(p(2),:) !< V1 = XVirt-X2
			V1(3) = 1.0d0

			V2(1:2) = Coordinates(p(3),:)-Coordinates(p(2),:) !< V2 = X3-X2
			V2(3) = 0.0d0

			!< This is already a normal vector
			V1 = Cross(V1, V2) !< Cross vector product is not assumed to have the vectors V1 and V2 well positioned

			Xm = Xm - Coordinates(p(1),:) !< Outer vector of edge oposed to vertex 1
			if (dot_product(Xm,V1(1:2))>0) then
				this%OuterVector(p(1),:) = V1(1:2)
			else
				this%OuterVector(p(1),:) = -V1(1:2)
			end if
			call Normalize(this%OuterVector(p(1),:))
        end do
    end subroutine



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Common Use on the manipulation of TriangularElements!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


    !> \brief Base functions on Triangular Linear coordinate $( \eta_1, \eta_2, \eta_3 )$ system
    !!
    !! \param i - the index of the base function
    !! \param eta - double precision 3d vector of the point on which the base is calculated
    !! \return v - double precision scalar calculated value
    !!
    pure function TriangularPhi(i, eta ) result (v)
        implicit none
        integer, intent(in) :: i
        real(kind=8), dimension(3), intent(in) :: eta
        real(kind=8) :: v
		if(size(eta) /= 3 ) return
        select case(i)
            case(1)
                v = eta(1)
            case(2)
                v = eta(2)
            case(3)
                v = eta(3)
            case default
                v = 0.0d0
        end select
    end function


    !> \brief Base functions on single variable, or linear, system
    !!
    !! \param i - the index of the base function
    !! \param eta - double precision scalar point on which the function is calculated
    !! \return v - double precision scalar calculated value
    !!
    pure function LinearPhi(i, eta) result (v)
        implicit none
        integer, intent(in) :: i
        real(kind=8), intent(in) :: eta
        real(kind=8) :: v

        select case(i)
            case(1)
                v = 0.5d0 * (1.0d0 - eta)
            case(2)
                v = 0.5d0 * (1.0d0 + eta)
            case default
                v = 0.0d0
        end select
    end function

    !> \brief Returns the Jacobian over the 1-D Boundary, or side of the Triangle, of integrals performed
    !!
    !! \param X - the Nodes Coordinates
    !! \param indexes - 3D integer vector comprising the indexex on the rotation of sides on which calculate the Jacobian
    !! \return v - double precision value of the Jacobian
    !!
    pure function JacobianLength(X, indexes) result (v)
		implicit none
        real(kind=8), dimension(nNodes,probDim), intent(in) :: X
        integer, dimension(3), intent(in) :: indexes
        real(kind=8) :: v
        v = 0.5d0*( (X(indexes(3),1) - X(indexes(2),1)**2.0d0 + (X(indexes(3),2) - X(indexes(2),2)**2.0d0)**0.5d0
    end function


    !> \brief Returns a analytically calculated Matrix which comprises of the Temperature contribution within the element area
    !!
    !! \param this - the Element structure
    !! \param X - the nodes coordinates
    !! \param orientation - the Axis representing the Allocation of the Matrix
    !! \return local - The Analytically calculated matrix
    !!
	pure function MatrixH_l( this, X, orientation ) result (local)
		implicit none
		type(TriangleElement), intent(in) :: this
		real(kind=8), dimension(nNodes,probDim), intent(in) :: X
		integer, intent(in) :: orientation

		real(kind=8), dimension(nNodes,nNodes) :: local
!// \todo (ronaldo#1#): É Necessario adequar para a condutividade térmica em função da temperatura.

		select case(orientation)
			case(XAxis)
				local (1,1) = X(2,2) - X(3,2) !< $y_{23}$
				local (1,2) = local(1,1)
				local (1,3) = local(1,1)
				local (2,1) = X(3,2) - X(1,2) !< $y_{31}$
				local (2,2) = local(2,1)
				local (2,3) = local(2,1)
				local (3,1) = X(1,2) - X(2,2) !< $y_{12}$
				local (3,2) = local(3,1)
				local (3,3) = local(3,1)
			case(YAxis)
				local (1,1) = X(3,1) - X(2,1) !< $x_{32}$
				local (1,2) = local(1,1)
				local (1,3) = local(1,1)
				local (2,1) = X(1,1) - X(3,1) !< $x_{13}$
				local (2,2) = local(2,1)
				local (2,3) = local(2,1)
				local (3,1) = X(2,1) - X(1,1) !< $x_{21}$
				local (3,2) = local(3,1)
				local (3,3) = local(3,1)
			case default
				local (1,1) = X(2,2) - X(3,2) !< $y_{23}$
				local (1,2) = local(1,1)
				local (1,3) = local(1,1)
				local (2,1) = X(3,2) - X(1,2) !< $y_{31}$
				local (2,2) = local(2,1)
				local (2,3) = local(2,1)
				local (3,1) = X(1,2) - X(2,2) !< $y_{12}$
				local (3,2) = local(3,1)
				local (3,3) = local(3,1)

		end select
		local = local/6.0d0
	end function


    !> \brief Generate a element of a Matrix $E_{lr,i}$, or a Matrix $E_{lr,i,B}$
    !!  Element boundary thermal conductivity contribution on side $i$
    !!  To gerenate a $E_{lr,i,B}$ matrix usa a negative Stabilization parameter
    !! \param this -> the element
    !! \param l, r -> ranges from x to y, each
    !! \param i -> the boundary
    !! \param k, m -> Matrix internal indexes
    !! \param ThermCond -> Thermal Conductivity. Function which defines the Thermal Conductivity
    !! \param Stabilization -> Stabilization Parameter
    !! \return cell - double precision scalar
    !!
	pure function MatrixE_lr_i_km( this, l, r, i, k, m, X, ThermCond ) result (cell)
        implicit none
		integer, parameter :: nInterpolationPoints = 3
		integer, parameter :: nDimension = 2
        integer, parameter :: IntPoints = 1
        integer, parameter :: nSides = 3
        type(TriangleElement), intent(in) :: this
		integer, intent(in) :: l, r, i, k, m
        real(kind=8), dimension(nNodes,2), intent(in) :: X !< Planar Element
        interface
			pure function ThermCond(T) result (v) !< Thermal Conductivity Function
				implicit none
				real(kind=8), intent(in) :: T
				real(kind=8) :: v
			end function
		end interface

        real(kind=8) :: cell
		integer :: g !< integration points index
		integer :: a !< Thermal Conductivity interpolation Index

		real(kind=8), dimension(IntPoints) :: IntegrationPoints
		real(kind=8), dimension(IntPoints) :: IntegrationWeights

		integer, dimension(3) :: vertex

		IntegrationPoints = LinearIntegrationPoints(IntPoints)
		IntegrationWeights = LinearIntegrationWeights(IntPoints)

		cell = 0.0d0
		vertex = rotation(i)
		forall (a=2:3) !< Points not regarding the vertex pivoting the situation
			forall (g=1:IntPoints) !< Warp through Integration Points
				cell = cell + ThermCond(this%DOF(Temperature_Initial+a-1))*&
							LinearPhi(vertex(a), IntegrationPoints(g)) *&
							LinearPhi(k, IntegrationPoints(g)) *&
							LinearPhi(m, IntegrationPoints(g)) *&
							IntegrationWeights(g)
			end forall
		end forall
		cell = this%OuterVector(i,l) * this%OuterVector(i,r) * JacobianLength(X, vertex) * Stabilization * cell
	end function


    !> \brief Returns a Component cell of the Matrix $H_{l,i}$ Calculated over the side $i$ of the element
    !!
    !! \param this - the Element Structure
    !! \param l - the orientation Axis of the Calculated Matrix
    !! \param i - the side on which this calculations are being performed
    !! \param k - the local matrix $H_{l,i}$ line of the component
    !! \param m - the local matrix $H_{l,i}$ column of the component
    !! \param X - the Nodes components
    !! \param ThermCond - the function defining the Thermal Conductivity at node $n$ depending, only, on temperature $T_n$
    !! \return cell - double precision scalar
    !!
	pure function MatrixH_l_i_km( this, l, i, k, m, X, ThermCond ) result (cell)
        implicit none
		integer, parameter :: nInterpolationPoints = nNodes
		integer, parameter :: nDimension = probDim
        integer, parameter :: IntPoints = 1
        type(TriangleElement), intent(in) :: this
		integer, intent(in) :: l, i, k, m
        real(kind=8), dimension(nNodes,2), intent(in) :: X !< Planar Element
        interface
			pure function ThermCond(T) result (v) !< Thermal Conductivity Function
				implicit none
				real(kind=8), intent(in) :: T
				real(kind=8) :: v
			end function
		end interface

        real(kind=8) :: cell
		integer :: g !< integration points index
		integer :: a !< Thermal Conductivity interpolation Index

		real(kind=8), dimension(IntPoints,3) :: IntegrationPoints
		real(kind=8), dimension(IntPoints) :: IntegrationWeights

		IntegrationPoints = LinearIntegrationPoints(IntPoints)
		IntegrationWeights = LinearIntegrationWeights(IntPoints)

		cell = 0.0d0
		vertex = rotation(i)
		forall (a=2:3) !< Points not regarding the vertex pivoting the situation
			forall (g=1:IntPoints) !< Warp through Integration Points
				cell = cell + ThermCond(this%DOF(Temperature_Initial+a-1))*&
							LinearPhi(vertex(a), IntegrationPoints(g)) *&
							LinearPhi(k, IntegrationPoints(g)) *&
							LinearPhi(m, IntegrationPoints(g)) *&
							IntegrationWeights(g)
			end forall
		end forall
		cell = cell  * JacobianLength(X, vertex) * this%OuterVector(i,l)
	end function


    !> \brief Returns a component of the Matrix $J_{l,i}$ calculated over the side $i$ of the element
    !!
    !! \param this - the Element Structure
    !! \param l - the orientation Axis of the Calculated Matrix
    !! \param i - the side on which this calculations are being performed
    !! \param k - the local matrix $J_{l,i}$ line of the component
    !! \param m - the local matrix $J_{l,i}$ column of the component
    !! \param X - the Nodes components
    !! \return cell - double precision scalar
    !!
	function MatrixJ_l_i_km( this, l, i, k, m, X ) result (cell)
        implicit none
		integer, parameter :: nInterpolationPoints = 3
		integer, parameter :: nDimension = 2
        integer, parameter :: IntPoints = 1
        type(TriangleElement), intent(inout) :: this
		integer, intent(in) :: l, i, k, m
        real(kind=8), dimension(nNodes,2), intent(in) :: X !< Planar Element

        real(kind=8) :: cell
		integer :: g !< integration points index

		real(kind=8), dimension(IntPoints) :: IntegrationPoints
		real(kind=8), dimension(IntPoints) :: IntegrationWeights

		integer, dimension(3) :: vertex

		IntegrationPoints = LinearIntegrationPoints(IntPoints)
		IntegrationWeights = LinearIntegrationWeights(IntPoints)

		cell = 0.0d0
		vertex = rotation(i)
		forall (g=1:IntPoints) !< Warp through Integration Points
			cell = cell + LinearPhi(k, IntegrationPoints(g)) *&
						LinearPhi(m, IntegrationPoints(g)) *&
						IntegrationWeights(g)
		end forall
		cell = cell * JacobianLength(X, vertex) * this%OuterVector(i,l)
	end function


    !> \brief Analitically calculated Matrix component of the block-matrices $G_{T,i}$ and $G_{T,i,B}$
    !!
    !! \param i - the side on which the calculation is being performed
    !! \param X - The Nodes Coordinates
    !! \return C - The Matrix itself
    !!
   pure function MatrixC_i( i, X ) result (C)
        implicit none
        integer, intent(in) :: i
		real(kind=8), dimension(nNodes,probDim), intent(in) :: X

        real(kind=8), dimension(nNodes,nNodes) :: C

        integer, dimension(nNodes) :: vertex
        integer :: a, b

        C = 0.0d0
        vertex = rotation(i)
		C(vertex(2),vertex(2)) = 2.0d0
		C(vertex(3),vertex(3)) = 2.0d0
		C(vertex(2),vertex(3)) = 1.0d0
		C(vertex(3),vertex(2)) = 1.0d0
        C = C * ( &
					(X(vertex(2),1)-X(vertex(3),1))**2.0d0 + &
					(X(vertex(2),2)-X(vertex(3),2))**2.0d0 &
					)**0.5
    end function




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Use on the calculation of the Stiffness Matrix !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !> \brief Auxiliary function S, which gives the ponderation value for the calculations of the fixed-form, or analytically calculated matrices on the triangle
    !!
    !! \param X - the Nodes coordinates
    !! \return local - the value
    !!
	pure function S (X) result (local)
		implicit none
		real(kind=8), dimension(nNodes, probDim), intent(in) :: X
		real :: local
		local = ( (X(1,1)-X(3,1))*(X(2,2)-X(3,2)) - (X(2,1)-X(3,1))*(X(1,2)-X(3,2)) )/2.0d0
	end function


    !> \brief Returns the Elementary matrix E, refered to the block-diagonal of the Stiffness local Matrix
    !!
    !! \param X - the Nodes coordinates
    !! \return local - The Entire Matrix
    !!
	pure function Matrix_E(X) result (local)
		implicit none
		real(kind=8), dimension(nNodes,probDim), intent(in) :: X
		real(kind=8), dimension(nNodes,nNodes) :: local
		local = 1.0d0
		local(1,1) = 2.0d0; local(2,2) = 2.0d0; local(3,3) = 2.0d0
		local = local * S(X)/12.0d0
	end function


    !> \brief Set the E central matrices of the local matrix
    !! Matrix of the inner element contribution relative to the heat flux
    !! The E matrices are dependent only of a determinant obtained from the Vertices Coordinates
    !! The E matrices compond the diagonal blocks of the first assembly
    !! \param this - The Triangle Element itself
    !! \param S - The determinant parameter
    !!
    subroutine SetMatricesE( Matrix, X )
        implicit none
		real(kind=8), dimension(MatrixSize,MatrixSize), intent(out) :: Matrix
		real(kind=8), dimension(nNodes,probDim), intent(in) :: X

		integer :: l, i
		integer, dimension(nNodes) :: idx

		forall(l = 1:probDim, i=i:nNodes)
			forall(i=1:nNodes) !< Set Range for Matrices E
				idx(i) = nNodes*(l-1)+i+1
			end forall
			idx = SetRange(l)
			Matrix(idx,idx) = Matrix(idx,idx) + Matrix_E(X)
		end forall
    end subroutine


    !> \brief Set the $H_k$ matrix of the local matrix, $k = x,y$
    !! Matrix of the inner element contribution relative to Temperature
    !! \param this - The Triangular Element itself
    !! \param X - The Coordinates of the nodes
    !!
	subroutine SetMatricesH_l( this, Matrix, X )
		implicit none
		type(TriangleElement), intent(inout) :: this
		real(kind=8), dimension(MatrixSize, MatrixSize), intent(out) :: Matrix
		real(kind=8), dimension(nNodes,probDim), intent(in) :: X
		integer :: l

		forall( l = 1:probDim)
			Matrix(l:(l*nNodes),Temperature_Initial:Temperature_Final) = Matrix(l:(l*nNodes),Temperature_Initial:Temperature_Final) + MatrixH_l( this, X, l)
		end forall

	end subroutine


    !> \brief Set the $J_k$ matrix of the local matrix, $k = x,y$
    !! Matrix of the inner element contribution relative to Temperature
    !! \param this - The Triangular Element itself
    !! \param X - The Coordinates of the nodes
    !!
	subroutine SetMatricesJ_l( this, Matrix, X )
		implicit none
		type(TriangleElement), intent(inout) :: this
		real(kind=8), dimension(MatrixSize, MatrixSize), intent(out) :: Matrix
		real(kind=8), dimension(nNodes,probDim), intent(in) :: X
		integer :: l

		forall( l = 1:probDim)
			Matrix(l:(l*nNodes),Temperature_Initial:Temperature_Final) = Matrix(l:(l*nNodes),Temperature_Initial:Temperature_Final) + MatrixH_l( this, X, l)
		end forall

	end subroutine


	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Face components of the Stiffness Matrix!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    !> \brief Set the $E_{lr,i}$ matrices comprising the element contribuitions to the Stiffness matrix on every side of the element
    !!
    !! \param this - the element on which the calculations will be performed
    !! \param Matrix - the Matrix on which the calculations will be recorded. It can be the element local matrix or a equivalent matrix for a external use
    !! \param X - the nodes coordinates
    !! \param ThermCond - The Thermal Conductivity function which determines the isotropic thermal conductivity coeficient
    !! \param Stabilization - double precision scalar for stabilization uses on the calculations
    !!
    subroutine SetMatricesE_lr_i( this, Matrix, X, ThermCond, Stabilization )
        implicit none
		integer, parameter :: nInterpolationPoints = nNodes
		integer, parameter :: nDimension = probDim
        integer, parameter :: nSides = 3
        type(TriangleElement), intent(inout) :: this
        real(kind=8), dimension(MatrixSize,MatrixSize+1), intent(inout) :: Matrix
        real(kind=8), dimension(nNodes,2), intent(in) :: X !< Planar Element
        interface
			pure function ThermCond(T) result (v) !< Thermal Conductivity Function
				implicit none
				real(kind=8), intent(in) :: T
				real(kind=8) :: v
			end function
		end interface
		real(kind=8), intent(in) :: Stabilization

		integer :: l, r  !< l, r = x, y, z (1, 2, 3)
        integer :: i !< side index
        integer :: k, m !< local base index
		integer :: l_m, r_m  !< local matrix real position

        forall(i = 1:nSides)  !< Warp through  Triangle1's Sides
			forall (l=1:nDimension) !< Warp through Local dimensions
				forall (r=1:nDimension) !< Warp through Local dimensions
					forall (k=1:nInterpolationPoints) !< Warp through  local basis
						forall (m=1:nInterpolationPoints) !< Warp throug local basis
							l_m = (l-1)*nNodes+k  !< Matrix Index of Element Matrix
							r_m = (r-1)*nNodes+m  !< Matrix Index of Element Matrix
							Matrix(l_m, r_m) = Matrix(l_m, r_m) + &
									MatrixE_lr_i_km(this, l, r, i, k, m, X, ThermCond, Stabilization) *&
									Stabilization
						end forall
					end forall
				end forall
			end forall
        end forall
    end subroutine


    !> \brief Sets the $H_{l,i}$ matrices comprising the Contributions of Temperature on each side of the Element
    !!
    !! \param this - the element on which the calculations will be performed
    !! \param Matrix - the Matrix on which the calculations will be recorded. It can be the element local matrix or a equivalent matrix for a external use
    !! \param X - the nodes coordinates
    !! \param ThermCond - The Thermal Conductivity function which determines the isotropic thermal conductivity coeficient
    !! \param Stabilization - double precision scalar for stabilization uses on the calculations
    !!
	subroutine SetMatricesH_l_i( this, Matrix, X, ThermCond, Stabilization )
        implicit none
		integer, parameter :: nInterpolationPoints = nNodes
		integer, parameter :: nDimension = probDim
        integer, parameter :: nSides = 3
        type(TriangleElement), intent(inout) :: this
        real(kind=8), dimension(MatrixSize, MatrixSize+1), intent(inout) :: Matrix
        real(kind=8), dimension(nNodes,2), intent(in) :: X !< Planar Element
        interface
			pure function ThermCond(T) result (v) !< Thermal Conductivity Function
				implicit none
				real(kind=8), intent(in) :: T
				real(kind=8) :: v
			end function
		end interface
        real(kind=8), intent(in) :: Stabilization

		integer :: l_m, r_m
		integer :: i
		integer :: l, r
		integer :: k, m

		forall(i = 1:nSides)
			forall(l=1:nDimension)
				forall(r=1:nDimension)
					forall(k=1:nInterpolationPoints)
						forall(m=1:nInterpolationPoints)
							l_m = (l-1)*nNodes+k
							r_m = (r-1)*nNodes+m
							Matrix(l_m, r_m) = Matrix(l_m, r_m) + ( -(0.5d0 + Stabilization)) * &
									MatrixH_l_i_km(this, l, i, k, m, X, ThermCond)
						end forall
					end forall
				end forall
			end forall
		end forall
	end subroutine


    !> \brief Sets the $J_{l,i}$ Temperatures comprising the Heat Fluxed contribution to the temperature on each side of the element
    !!
    !! \param this - the element on which the calculations will be performed
    !! \param Matrix - the Matrix on which the calculations will be recorded. It can be the element local matrix or a equivalent matrix for a external use
    !! \param X - the nodes coordinates
    !! \param Stabilization - double precision scalar for stabilization uses on the calculations
    !!
	subroutine SetMatricesJ_l_i( this, Matrix, X, Stabilization )
        implicit none
		integer, parameter :: nInterpolationPoints = nNodes
		integer, parameter :: nDimension = probDim
        integer, parameter :: nSides = 3
        type(TriangleElement), intent(inout) :: this
        real(kind=8), dimension(MatrixSize, MatrixSize+1), intent(inout) :: Matrix
        real(kind=8), dimension(nNodes,2), intent(in) :: X !< Planar Element
        real(kind=8), intent(in) :: Stabilization

		integer :: l_m, r_m
		integer :: i
		integer :: l, r
		integer :: k, m

		forall(i = 1:nSides)
			forall(l=1:nDimension)
				forall(r=1:nDimension)
					forall(k=1:nInterpolationPoints)
						forall(m=1:nInterpolationPoints)
							l_m = (l-1)*nNodes+k
							r_m = (r-1)*nNodes+m
							Matrix(l_m, r_m) = Matrix(l_m, r_m) + ( -(0.5d0 - Stabilization)) * &
									MatrixJ_l_i_km(this, l, i, k, m, X, Stabilization)
						end forall
					end forall
				end forall
			end forall
		end forall
	end subroutine


    !> \brief Sets the $G_{T,i}$ matrices of the Stiffness Matrix comprising the Temperature Contributions to the Temperature itself on each side
    !!
    !! \param Matrix - The matrix which the calculations will be recorded
    !! \param X - the Nodes coordinates
    !! \param Stabilization - double precision scalar Stabilization parameter
    !!
    subroutine SetMatrixG_T_i( Matrix, X, Stabilization )
        implicit none
        integer, parameter :: nSides = nNodes
		real(kind=8), dimension(MatrixSize,MatrixSize), intent(out) :: Matrix
		real(kind=8), dimension(nNodes,probDim), intent(in) :: X

		integer :: l

		forall(l = 1:nSides)
			Matrix(idx,idx) = Matrix + Stabilization * MatrixC_i(l, X)
		end forall
    end subroutine


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!! Use on the calculations of the Neighbour dominated vector of the Element !!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



    !> \brief Set the $E_lr_i_B$ comprising the Neighbours contribution to the independent constant vector due to the temperature effect on the flux
    !!
    !! \param this - the Element on which the computation is being performed
    !! \param Vector - the Vector on which the computations will be stored
    !! \param Neighbours - the Element Neighbours whom will give the Stabilized, or Committed, DoF
    !! \param X - the Nodes coordinates
    !! \param Stabilization - double precision scalar Stabilization parameter
    !! \param ThermCond - The Thermal Conductivity function for update the values
    !!
    subroutine SetMatricesE_lr_B_i( this, Vector, Neighbours, X, ThermCond, Stabilization )
        implicit none
		integer, parameter :: nInterpolationPoints = nNodes
		integer, parameter :: nDimension = probDim
        integer, parameter :: nSides = 3
        type(TriangleElement), intent(inout) :: this
        real(kind=8), dimension(MatrixSize), intent(inout) :: Vector
        type(TriangleElement), dimension(nSides), intent(in) :: Neighbours
        real(kind=8), dimension(nNodes,2), intent(in) :: X !< Planar Element
        interface
			pure function ThermCond(T) result (v) !< Thermal Conductivity Function
				implicit none
				real(kind=8), intent(in) :: T
				real(kind=8) :: v
			end function
		end interface
		real(kind=8), intent(in) :: Stabilization

		integer :: l, r  !< l, r = x, y, z (1, 2, 3)
        integer :: i !< side index
        integer :: k, m !< local base index
		integer :: l_m, r_m  !< local matrix real position

        forall(i = 1:nSides)  !< Warp through  Triangle1's Sides
			forall (l=1:nDimension) !< Warp through Local dimensions
				forall (r=1:nDimension) !< Warp through Local dimensions
					forall (k=1:nInterpolationPoints) !< Warp through  local basis
						forall (m=1:nInterpolationPoints) !< Warp throug local basis
							l_m = (l-1)*nNodes+k  !< Matrix Index of Element Matrix
							r_m = (r-1)*nNodes+m  !< Matrix Index of Element Matrix
							Vector(l_m) = Vector(l_m) - &
									MatrixE_lr_i_km(this, l, r, i, k, m, X, ThermCond, Stabilization) *&
									Stabilization * Neighbours(i)%DoF(r_m, Dof_Commit) * (-1.0d0) ! -1 is due the change on the equality sign
						end forall
					end forall
				end forall
			end forall
        end forall
    end subroutine


    !> \brief Sets the matrices $H_{l,B,i}$ comprising the Neighbours contribution to the independent constant vector due to the temperature effect on the flux
    !!
    !! \param this - the Element on which the computation is being performed
    !! \param Vector - the Vector on which the computations will be stored
    !! \param Neighbours - the Element Neighbours whom will give the Stabilized, or Committed, DoF
    !! \param X - the Nodes coordinates
    !! \param Stabilization - double precision scalar Stabilization parameter
    !! \param ThermCond - The Thermal Conductivity function for update the values
    !!
	subroutine SetMatricesH_l_B_i( this, Vector, Neighbours, X, Stabilization, ThermCond )
        implicit none
		integer, parameter :: nInterpolationPoints = nNodes
		integer, parameter :: nDimension = probDim
        integer, parameter :: nSides = 3
        type(TriangleElement), intent(inout) :: this
        real(kind=8), dimension(MatrixSize), intent(inout) :: Vector
        type(TriangleElement), dimension(nSides), intent(in) :: Neighbours
        real(kind=8), dimension(nNodes,2), intent(in) :: X !< Planar Element
        real(kind=8), intent(in) :: Stabilization
        interface
			pure function ThermCond(T) result (v) !< Thermal Conductivity Function
				implicit none
				real(kind=8), intent(in) :: T
				real(kind=8) :: v
			end function
		end interface

		integer :: l_m, r_m
		integer :: i
		integer :: l, r
		integer :: k, m

		forall(i = 1:nSides)
			forall(l=1:nDimension)
				forall(r=1:nDimension)
					forall(k=1:nInterpolationPoints)
						forall(m=1:nInterpolationPoints)
							l_m = (l-1)*nNodes+k
							r_m = (r-1)*nNodes+m
							Vector(l_m) = Matrix(l_m, r_m) + ( -(0.5d0 - Stabilization)) * &
									MatrixH_l_i_km(this, l, i, k, m, X, ThermCond) * &
									Neighbours(i)%DoF(r_m, Dof_Committ) * (-1.0d0) ! -1 is due the change on the equality sign
						end forall
					end forall
				end forall
			end forall
		end forall
	end subroutine


    !> \brief Sets the matrices $J_{l,B,i}$ comprising the Neighbours contribution to the independent constant vector due to the flux effect on the temperature
    !!
    !! \param this - the Element on which the computation is being performed
    !! \param Vector - the Vector on which the computations will be stored
    !! \param Neighbours - the Element Neighbours whom will give the Stabilized, or Committed, DoF
    !! \param X - the Nodes coordinates
    !! \param Stabilization - double precision scalar Stabilization parameter
    !!
	subroutine SetMatricesJ_l_B_i( this, Vector, Neighbours, X, Stabilization )
        implicit none
		integer, parameter :: nInterpolationPoints = nNodes
		integer, parameter :: nDimension = probDim
        integer, parameter :: nSides = 3
        type(TriangleElement), intent(inout) :: this
        real(kind=8), dimension(MatrixSize), intent(inout) :: Vector
        type(TriangleElement), dimension(nSides), intent(in) :: Neighbours
        real(kind=8), dimension(nNodes,2), intent(in) :: X !< Planar Element
        real(kind=8), intent(in) :: Stabilization

		integer :: l_m, r_m
		integer :: i
		integer :: l, r
		integer :: k, m

		forall(i = 1:nSides)
			forall(l=1:nDimension)
				forall(r=1:nDimension)
					forall(k=1:nInterpolationPoints)
						forall(m=1:nInterpolationPoints)
							l_m = (l-1)*nNodes+k
							r_m = (r-1)*nNodes+m
							Vector(l_m) = Vector(l_m) + ( -(0.5d0 + Stabilization)) * &
									MatrixJ_l_i_km(this, l, i, k, m, X) * &
									Neighbours(i)%DoF(r_m, Dof_Committ) * (-1.0d0) ! -1 is due the change on the equality sign
						end forall
					end forall
				end forall
			end forall
		end forall
	end subroutine


    !> \brief Sets the matrices $G_{l,B,i}$ comprising the Neighbours contribution to the independent constant vector due to the flux effect on the temperature
    !!
    !! \param Vector - the Vector on which the computations will be stored
    !! \param Neighbours - the Element Neighbours whom will give the Stabilized, or Committed, DoF
    !! \param X - the Nodes coordinates
    !! \param Stabilization - double precision scalar Stabilization parameter
    !!
    subroutine SetMatrixG_T_B_i( Vector, Neighbours, X, Stabilization )
        implicit none
        real(kind=8), dimension(MatrixSize), intent(inout) :: Vector
        type(TriangleElement), dimension(nSides), intent(in) :: Neighbours
		real(kind=8), dimension(nNodes,probDim), intent(in) :: X

		integer :: l

		forall(l = 1:nSides)
			Matrix(idx,idx) = Matrix - Stabilization * MatrixC_i(l, X) * &
									Neighbours(i)%DoF(r_m, Dof_Committ) * (-1.0d0) ! -1 is due the change on the equality sign
		end forall
    end subroutine




end module
